function [B, B_deriv]=construct_splines(x_interest, xi, n, k, p_low, p_high, down)

N=size(x_interest,1);

if (k>1),

    B=NaN*ones(N,p_high-p_low+1);
    B_deriv=B;
    [B_prev, B_prev_deriv]=construct_splines(x_interest, xi, n, k-1,-(k+down),n+down,down+1);
    
    for pp=p_low:p_high,
        pp_mod=pp-p_low+1;     
        B(:,pp_mod) = ((x_interest - ones(N,1)*xi(pp_mod,1)).*B_prev(:,pp_mod) + (xi(pp_mod+k+1,1)*ones(N,1) - x_interest).*B_prev(:,pp_mod+1)) / (xi(pp_mod+k+1,1)-xi(pp_mod,1));
        B_deriv(:,pp_mod) = ((x_interest - ones(N,1)*xi(pp_mod,1)).*B_prev_deriv(:,pp_mod) + (xi(pp_mod+k+1,1)*ones(N,1) - x_interest).*B_prev_deriv(:,pp_mod+1) + B_prev(:,pp_mod) - B_prev(:,pp_mod+1)) / (xi(pp_mod+k+1,1)-xi(pp_mod,1));
    end;

end;

if (k==1),

    B=0*ones(N,p_high-p_low+1);
    B_deriv=B;
    
    for ii=p_low:p_high,
        ii_mod=ii - p_low+2;
        B(:,ii_mod-1)=B(:,ii_mod-1) + ((x_interest>=xi(ii_mod,1)).*(x_interest<xi(ii_mod+1,1))) .* (ones(N,1)*xi(ii_mod+1,1)-x_interest) / ((xi(ii_mod+1,1)-xi(ii_mod-1,1))*((xi(ii_mod+1,1)-xi(ii_mod,1))));
        B_deriv(:,ii_mod-1)=B_deriv(:,ii_mod-1) + ((x_interest>=xi(ii_mod,1)).*(x_interest<xi(ii_mod+1,1))) .* (-1) / ((xi(ii_mod+1,1)-xi(ii_mod-1,1))*((xi(ii_mod+1,1)-xi(ii_mod,1))));
    end;    

    for ii=p_low:p_high,
        ii_mod=ii - p_low+1;
        B(:,ii_mod)=B(:,ii_mod) + ((x_interest>=xi(ii_mod,1)).*(x_interest<xi(ii_mod+1,1))) .* (x_interest-ones(N,1)*xi(ii_mod,1)) / ((xi(ii_mod+1,1)-xi(ii_mod,1))*((xi(ii_mod+2,1)-xi(ii_mod,1))));       
        B_deriv(:,ii_mod)=B_deriv(:,ii_mod) + ((x_interest>=xi(ii_mod,1)).*(x_interest<xi(ii_mod+1,1))) .* (1) / ((xi(ii_mod+1,1)-xi(ii_mod,1))*((xi(ii_mod+2,1)-xi(ii_mod,1))));       
    end;    
        
end;

    %
